function [phi, dphidxi, dphideta, dphiddxi, dphidxideta, dphiddeta] = fe_basis(FE_Type, FE_Order, qx, qy)
% function [phi, dphidxi, dphideta, dphiddxi, dphidxideta, dphiddeta] =
% fe_basis(FE_Type, FE_Order, qp(1,:),qp(2,:));

one  =  ones(size(qx,1),1);
zero = zeros(size(qx,1),1);
r = qx;  s = qy; t = 1 - r - s;

% if only low order FE is used, all FEs are the same.
if(FE_Order == 1)
    phi         = [r,          s,          t];
    dphidxi   = [one , zero,   - one];
    dphideta = [zero,  one,   - one];
    dphiddxi = [zero, zero,  zero];
    dphidxideta = [zero, zero,  zero];
    dphiddeta = [zero, zero,  zero];
    return;
end


%%  type 1
if (strcmp(FE_Type, 'Lagrange') == 1)
    switch (FE_order)
        case 2
            phi      = [r.*(2*r-1),  s.*(2*s-1),  t.*(2*t-1),  4*s.*t,  4*t.*r,  4*r.*s];
            dphidxi  = [4*r-1,  zero, 1-4*t, -4*s, 4*(t-r), 4*s];
            dphideta = [zero, 4*s-1, 1-4*t, 4*(t-s), -4*r, 4*r];
            
        case 3
            
        case 4
            
        case 5
            
        otherwise
            error('Too high order FE for Lagrange Polynomial are not recommended!\n');
    end
            
%%  type 2
elseif(strcmp(FE_Type,  'BB') == 1)
    L1 = r;
    L2 = s;
    L3 = 1 - r - s;
    
    phi = vdm23_bb(FE_Order, L1, L2, L3);
    
    desc = desc_pattern(FE_Order);

    % first derivatives on master elements
    Id = diag(ones((FE_Order + 1)*(FE_Order + 2)/2, 1));
    u = FE_Order*de_cast_step(Id, FE_Order, 1, 0, -1, desc);  % direction derivetives
    v = FE_Order*de_cast_step(Id, FE_Order, 0, 1, -1, desc);

    Mat = vdm23_bb(FE_Order - 1, L1, L2, L3);  % low order vdm
    dphidxi  = Mat*u;  %get the deriv B-form value
    dphideta = Mat*v;

    % second derivatives on master elements
    uu = (FE_Order - 1)*de_cast_step(u, FE_Order - 1, 1, 0, -1, desc);
    uv = (FE_Order - 1)*de_cast_step(u, FE_Order - 1, 0, 1, -1, desc);
    vv = (FE_Order - 1)*de_cast_step(v,  FE_Order - 1, 0, 1, -1, desc);

    Mat = vdm23_bb(FE_Order - 2, L1, L2, L3);  % low order vdm
    dphiddxi     = Mat*uu;  %get the deriv B-form value
    dphidxideta = Mat*uv;
    dphiddeta    = Mat*vv;

% The following order is following FE
%     switch (FE_order)
%         case 2
%             psi      = [r.*r,  s.*s, t.*t,  2*s.*t,  2*t.*r,  2*r.*s];
%             psi_xi  = [2*r,  zero, -2*t, -2*s, 2*(t-r), 2*s];
%             psi_eta = [zero, 2*s, -2*t, 2*(t-s), -2*r, 2*r];
%             
%         otherwise
%             
%     end

    
%%  only for FangISO2, just for test
elseif(strcmp(FE_Type,  'Hierarchy') == 1)
    switch (FE_order)
        case 2
            phi          = [r,  s,  t,             4*s.*t,  4*r.*t,  4*r.*s];
            dphidxi   = [one,zero,-one,     -4*s, 4*(t-r),      4*s];
            dphideta  = [zero,one,-one, 4*(t-s),    -4*r,      4*r];
            
        otherwise
            error('Sorry, the FE type is not supported!\n');
    end
    
%     %%     type 3   
%     elseif(strcmp(FE_Type, 'DG') == 1)
%         switch (FE_order)
%             case 2   % assume the v1(0,0). v2(1,0), v3(0,1) for the  master element
%                 psi = [r.*(2*r - 1), s.*(2*s - 1), t.*(2*t-1), 4*s.*t, 4*t.*r, 4*r.*s];
%                 psi_xi  = [4*r - 1, zero, -3 + 4*r + 4*s, -4*s,  4 - 8*r - 4*s, 4*s];
%                 psi_eta = [zero, 4*s-1, -3 + 4*s + 4*r, 4 - 8*s - 4*r, -4*r, 4*r];
% 
%             otherwise
%         end
% 
%%  do not supported warning
else
    error('Sorry, the FE type is not supported!\n');
end

end

%%%%%%%%%
%
%
function Mat = vdm23_bb(d,b1,b2,b3)
% comput the Bform of degreed d at points b1,b2,b3;
m = (d+1)*(d+2)/2;
plot_m = length(b1);
[I,J,K] = indices(d);
IM = diag(I)*ones(m,plot_m);
JM = diag(J)*ones(m,plot_m);
KM = diag(K)*ones(m,plot_m);

plot_IM = diag(b1)*ones(plot_m,m);
plot_JM = diag(b2)*ones(plot_m,m);
plot_KM = diag(b3)*ones(plot_m,m);
Mat = (plot_IM).^(IM').*(plot_JM).^(JM').*(plot_KM).^(KM');
IF = gamma(I+1);
JF = gamma(J+1);
KF = gamma(K+1);
A = factorial(d)*ones(plot_m,m)*diag(1./(IF.*JF.*KF));
Mat = A.*Mat;
end

%%%%%%%%%
%
%
function desc = desc_pattern(d)
% this function may be execute only once in whole comput session.
% much like the function asce_pattern.
m = (d+1)*(d+2)/2;
desc = zeros(m,3);
begin = 1;
for j = 0:d
    idx = (begin:(begin+j))';
    desc(idx,1) = idx;
    desc(idx,2) = idx + j + 1;
    desc(idx,3) = idx + j + 2;
    begin = begin + j + 1;
end
end

%%%%%%%%%
%
%
function Bout = de_cast_step(Bin,d,lam1,lam2,lam3,desc_pattern)
m_in = size(Bin,1);
m_out = m_in-d-1;
n = length(lam1);
indx1 = desc_pattern(1:m_out,1);  % always 1:m_out, but we place it here for simplicity
indx2 = desc_pattern(1:m_out,2);
indx3 = desc_pattern(1:m_out,3);
if size(Bin,2) == 1
   Bout = Bin(indx1)*lam1' + Bin(indx2)*lam2' + Bin(indx3)*lam3';
else
   Bout = Bin(indx1,:)*spdiags(lam1,0,n,n) + Bin(indx2,:)*spdiags(lam2,0,n,n) + ...
      Bin(indx3,:)*spdiags(lam3,0,n,n);
end
end